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1 Introduction 

The main disadvantages of Magnetic Resonance Imaging (MRI) are its long 
scan times and, in consequence, its sensitivity to motion. Exploiting the com¬ 
plementary information from multiple receive coils, parallel imaging is able to 
recover images from under-sampled k-space data and to accelerate the measure¬ 
ment 

mu Hum® 13 . Because parallel magnetic resonance imaging can be 
used to accelerate basically any imaging sequence it has many important appli¬ 
cations. Parallel imaging brought a fundamental shift in image reconstruction: 
Image reconstruction changed from a simple direct Fourier transform to the so¬ 
lution of an ill-conditioned inverse problem. This work gives an overview of im¬ 
age reconstruction from the perspective of inverse problems. After introducing 
basic concepts such as regularization, discretization, and iterative reconstruc¬ 
tion, advanced topics are discussed including algorithms for auto-calibration, 
the connection to approximation theory, and the combination with compressed 
sensing. 


2 Parallel Imaging as Inverse Problem 

2.1 Forward Model 

The signal from a receive-coil array with N channels is given by [5]: 

yj(t) = f d fm{f)Cj(r)e~ 2m ^' r 1 < j < N (1) 

J n 

The complex-valued magnetization image m represents the state of the trans¬ 
verse magnetization of the excited spins in the held-of-view (FOV) C at 
the time of image acquisition. In MRI, typically a volume (d = 3) or a thin 
slice (d = 2) of proton spins is excited using a resonant radio-frequency pulse. 
The image (or volume) is modulated by the complex-valued sensitivities Cj of 
all receive coils (Fig. [l]). The k-space signals yj(t) are then given by the Fourier 
transform of the coil images Cjrri sampled at discrete time points f; along a 
given k-space trajectory k(t). Equation [l] neglects relaxation and off-resonance 
effects during the acquisition, which is possible if all samples at time points ti 
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Figure 1: All individual coil images mcj for a phased-array coil with 32 elements 
(channels). For each coil element j = 1,... ,32, the complex-valued magneti¬ 
zation image m is modulated by its unique receive sensitivity Cj. The phase is 
color coded. 
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Figure 2: Under-sampling in k-space causes aliasing artifacts in the image do¬ 
main. Images of a human brain were reconstructed with a direct Fourier trans¬ 
form from fully-sampled data (full) and four-fold under-sampled data (regular 
under-sampling in one dimension (1x4), two dimensions (2x2), and Poisson-disc 
sampling). 
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Figure 3: Parallel imaging is an inverse problem. The forward model is given 
by a composition of a physical model and a sampling operator F = P o G. The 
operator G maps magnetization images m £ L(fl, C) on a field-of-view (FOV) 
to the ideal k-space signals / £ L(C d ,C N ) from N channels. The sampling 
operator P maps ideal signals / to k-space samples fj(k(ti)) for all channels 
j = 1,..., N and at all sample locations k(ti) £ S. 


are acquired in a small window around the echo time (TE) after excitation. Re¬ 
laxation and off-resonance effects from the excitation of the spins until the echo 
time are incorporated into the image and define the image contrast. Because 
in most cases it is not possible to acquire all data in this short acquisition win¬ 
dow, samples have to be acquired in a repeated series of identical experiments, 
which always restore the magnetization image to exactly the same stateQ This 
requirement to repeat the basic experiment many times is the reason for the 
long scans times in MRI. The goal of parallel imaging is to reduce the amount 
of data required to reconstruct images by optimally exploiting the complemen¬ 
tary information from multiple receive coils. Although there are fundamental 
limits to the encoding power of the receive sensitivities, it has the potential to 
accelerate MRI by a factor of about four in each spatial dimension [5]. 

2.2 Image Reconstruction 

If the coil sensitivities are known (e.g. from a pre-scan), image reconstruction 
for parallel imaging can be formulated as a linear inverse problem with discrete 
data m ■ Mathematically, the forward problem is given by an operator F which 
maps the magnetization image m £ L(fl, C) of excited spins in a FOV fl c R d 
to the sample values yj(ti). This operator can be thought of as the composition 

1 Single-shot Echo-Planar Imaging (EPI) and spiral imaging sequences acquire all data in 
a single acquisition. These sequences are fast, but image quality is compromised by blur¬ 
ring, distortions, and phase cancellation, due to relaxation and off-resonance effects. Parallel 
imaging can be used to shorten the acquisition window and reduce these artifacts. 
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Figure 4: Images of a human brain reconstructed using parallel imaging from 
fully-sampled data (full) and under-sampled data (acceleration by 2x2, 3x3, 
and 4x4) acquired with a 32-channel coil. Zi-wavelet regularization was used to 
suppress noise in the reconstruction. 


F = P o G of a physical model G for hypothetical multi-dimensional continuous 
k-space signals and a sampling operator P (see Fig.|3|. The operator G is given 

by 

G : L{n, C) -A C°°(R n , C) 

to i-A / with fj{k) = (m,enc?’ k ), (2) 

with the encoding functions defined as enc’’ k (r) := Cj(r)e +2mk '- t ' l ' r and a scalar 
product defined on L(f2,C) (anti-linear in the second argument). The sampling 
operator P evaluates the ideal k-space signals fj at the sample locations k(ti) G 
S in k-space. It is assumed that the sampling process corresponds to the point- 
evaluation of ideal k-space functions, i.e. yj{ti) = The sample values 

fj(S) are corrupted by additive white Gaussian noise. Although the noise is 
typically correlated between receive channels, this correlation can be removed 
with a whitening step. 

A variational solution to the inverse problem can be defined as the minimizer 
of a functional, i.e. 

m a := argmin^. \\Fx — y ||| + aR(x) . (3) 

The functional is composed of a least-squares data fidelity term (which alterna¬ 
tively may also include weighting or use a robust norm Hum]) and an additional 
regularization term R. Discretized versions of this minimization problem are the 
basis of SMASH and SENSE parallel imaging methods [5l l6l 171 IT3l fT4] . For par¬ 
allel MRI, this formulation has two advantages: First, arbitrary Cartesian or 
non-Cartesian sampling schemes can be used na. Second, the regularization 
term can be used to introduce prior knowledge about the solution^] Figure [4] 

2 A third advantage - which conceptually goes beyond parallel imaging - is the possibility 
to extend the forward model to include further physical effects in model-based reconstruction, 
e.g. field maps |151 I16| , motion-induced phase maps El E], motion El SS, relaxation 
maps |21ll22l[23| . or diffusion models [22], etc. 
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shows images of a human brain recovered from under-sampled data by numerical 
optimization of Equation [3] 

2.3 Regularization 

Ill-conditioning causes noise amplification during image reconstruction, which 
initially limited the application of parallel imaging to only moderate accelera¬ 
tion. This limitation can be overcome by incorporating prior knowledge about 
the image using regularization methods [2S1H31 [2S] - In the simplest case, regu¬ 
larization may consist of a basic quadratic penalty in the framework of a linear 
reconstruction, or make use of much more sophisticated techniques which exploit 
the structure of images but demand a non-linear reconstruction. For a least- 
squares problem with quadratic regularization, i.e. R{x ) = || \[W[x — rco) ||| with 
a positive definite operator W, the solution of Eq. [3] is explicitly given by the 
formula 


m a = x 0 + (F h F + aW)- 1 F H (y — Fx 0 ) . (4) 

S -V-' 

In the limit a —> 0 this solution is called the best approximate and is given 
by the Moore-Penrose pseudo-inverse F'. It has a statistical interpretation - 
assuming white Gaussian noise - as the best unbiased estimate for the image. 
Regularization can be interpreted as prior knowledge and the optimizer as a 
maximum a posteriori (MAP) estimate of the image. Although regularization 
leads to a fundamental trade-off between bias and noise - which has to be chosen 
carefully for optimal image quality - it makes the use of higher acceleration 
possible. An optimal estimate in terms of mean squared error can only be 
obtained with regularization. 

For optimal results, the prior knowledge should include as much specific 
knowledge about the image as possible. For example, regularization can exploit 
smoothness in the time domain m, or exploit that changes relative to a fixed 
reference image Xq can be assumed to be small. The later is used successfully for 
real-time MRI [28] or dynamic contrast enhanced (DCE) MRI [29]. While I 2 - 
regularization is simple to implement and already a clear improvement compared 
to unregularized parallel imaging, much better results can be be obtained when 
using more advanced techniques such as Zi-wavelet regularization, i.e. R(x) = 
|| DWT x||i, total variation, or other edge-preserving penalties m eh m sa. 

2.4 Discretization 

Numerical reconstruction methods make use of discretization, i.e. the unknown 
image is expanded into a sum of basis functions: 

m(r) ^^aigi(r) (5) 

l 
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Figure 5: Reconstruction from 16-fold under-sampled data using low, medium, 
and high amount of I 2 and Zi-wavelet regularization. Reconstructed images and 
difference images are shown. There is a trade-off between bias (residual artifacts, 
blurring) and noise depending on the amount of regularization. 
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Although the choice of this basis has subtle implications for results and in¬ 
terpretation, this topic has not drawn much attention]^] Most image-domain 
formulations based on SENSE use a grid of Dirac pulses to represent the im¬ 
age, because multiplication with the coil sensitivities in the forward model is 
then simply a point-wise multiplication. In contrast, k-space methods such as 
SMASH use a finite Fourier basis. 

Discretization has a regularizing effect, i.e. the discretized problem might 
have better condition than the continuous problem. In parallel imaging, this 
effect can often be seen in the area outside of the sampled k-space region. Ex¬ 
trapolation to these area causes high noise amplification [52]. A discretization 
scheme which excludes these degrees of freedom will be less affected by noise. On 
the other hand, a small basis leads to discretization errors because the solution 
can not be represented accurately. Both problems can be avoided by using fine 
discretization with a large number of basis function and explicit regularization 
to control noise amplification 1551 . 

It should be noted that for parallel imaging the ideal continuous solution of 
Eq. [3] can usually be computed almost perfectly 135] • Coil sensitivities are very 
smooth and can be approximated with a small number of Fourier coefficients 
(on an over-sampled FOV). The forward operator can then be understood as a 
convolution of the Fourier series of the image with a short filter. Because the 
acquired k-space data consists of a finite number of samples, also only a finite 
number of low-order Fourier coefficients from the infinite number of coefficients 
in the Fourier series of the image actually appear in the result of this convolu¬ 
tion. For quadratic regularization, a minimum-norm solution is obtained when 
the infinite number of remaining higher-order coefficients are set to zero. An im¬ 
plementation of the forward operator requires an a-periodic convolution which 
can be implemented efficiently using a fast Fourier transform (FFT) algorithm. 
In practice this differs from a conventional SENSE implementation only by using 
zero padding and in the exact interpretation of the recovered coefficients. 

For non-quadratic regularization, discretization errors may also arise in the 
implementation the regularization terms. In general, oversampling can be used 
to reduce these errors. The combination of non-quadratic regularization and 
oversampling can also avoid artifacts caused by truncation of the signal in the 
Fourier domain (Gibbs ringing) [57] . Finally, an important aspect related to 
discretization is a common error called an “inverse crime” [38i: When testing a 
reconstruction algorithms with simulated data, computing this data using the 
same discretization scheme as used for the reconstruction can result in highly 
misleading results. One possibility to avoid this error is the use of analytic 
phantoms 155] , 

2.5 Numerical Optimization 

For regular sampling schemes and quadratic regularization, a solution can be 
computed directly with matrix inversion, because the equations decouple into 

‘^Discussions in terms of “ideal voxel functions” (or “target voxel shapes”) can be found in 
earlier works 0CE3]. 
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small systems [7J. Although very efficient, this approach is not very flexible. 
Matrix-free iterative methods can be used instead to efficiently compute the 
solution for arbitrary sampling schemes |14j . Matrix-free methods are build 
from procedural implementations of the matrix-vector products y H > F H y and 
x i—^ Fx (or x i—^ F h Fx). For Cartesian sampling, these operations can be 
implemented using point-wise multiplications and FFT algorithms. For non- 
Cartesian sampling, efficient non-uniform fast Fourier transform (nuFFT) al¬ 
gorithms have been developed to estimate samples at arbitrary k-space loca¬ 
tions [40; SU [423 . Even more efficient algorithms can be designed when consid¬ 
ering the combined operator F H F which appears in the gradient of the least- 
squares data fidelity. For example, the effect of sampling in the Fourier domain 
can be computed exactly as a convolution with a truncated point-spread func¬ 
tion with the use of two zero-padded FFTs [43] • Overlap-add and overlap-save 
convolution algorithms can be used to exploit the compact representation of the 
coil sensitivities in the Fourier domain [44] , 

For quadratic regularization, an efficient iterative algorithm is the conjugate 
gradient method applied to the normal equations El: 

(F h F + aW) x = F H y . (6) 

It should be noted that the use of a density compensation as known from the 
direct gridding algorithm is neither required nor recommended]^] 

For Zi-regularization, the simplest (and slowest) reconstruction algorithm is 
iterative soft-thresholding [451 : 

z n =x n + aF H (y - Fx n ) (7) 

x n = T~ 1 r] X (Tz n ) (8) 

The first equation is a gradient descent step and the second update uses soft- 
thresholding 77 in a transform basis T, e.g. a discrete wavelet transform. This 
scheme converges slowly, but can be accelerated with the addition of a ravine 
step as in FISTA j45], 47.. Especially when using multiple convex penalties R n , 
a very flexible approach is an extension of the Alternating Direction Method of 
Multipliers (ADMM) [45] @2 E2 that can solve optimization problems of the 
form 


L 

aTgmm x ^2 R n (B n x) . (9) 

n= 1 

This approach is very flexible and has many advantages from an implementation 
point of view, because it splits the optimization into independent sub-problems. 

4 In (non-iterative) gridding, the density compensation is a diagonal matrix which approx¬ 
imates the inverse of F H F. Combined with the adjoint F H it yields an approximation of the 
pseudo-inverse, i.e. F H D « F'. Including a density compensation into an iterative optimiza¬ 
tion method produces solutions different from the optimal least-squares solution m- I.e., 
naively using F H D instead of F H as is sometimes suggested to improve the condition yields 
a different optimization problem: argmin^ \\VT)(Fx — y) ||| + olR{x) 


Many different kinds of regularization terms can easily be integrated if respective 
proximal operators of the form 


prox^{y) := argmin^, ^\\x - y\\% + R{x) 


( 10 ) 


are available in a computationally efficient form. For example, the proximal 
operator for the data fidelity term is the ^-regularized least-square inverse 
which can be computed efficiently with the methods of conjugate gradients. 
The proximal operator for Zi-regularization can be evaluated simply using soft- 
thesholding. Efficient implementations of many advanced algorithms which 
make use of parallel programming can be found Berkeley Advanced Recon¬ 
struction Toolbox (BART) fol| . 

3 Auto-calibration 

To obtain optimal results in parallel MRI, accurate and up-to-date informa¬ 
tion about the sensitivities of all receive coils is required. While approximate 
coil sensitivities can be computed from the geometry of the receive coils using 
the Bios Savart law, exact sensitivities depend on the loading of the coils and 
need to be determined with high accuracy during the actual measurement. A 
pre-scan can provide accurate calibration information, but this requires that 
experimental conditions stay exactly the same for the duration of the whole 
examination. Because this is not always guaranteed, auto-calibration methods 
have been developed which perform calibration using a small amount of addi¬ 
tional data acquired during each individual scan (55, 55. [53J S3 ■ Because this 
reduces overall acceleration, optimal calibration from a minimum amount of 
data is desired. Two advanced techniques are described in the following: Joint 
estimation techniques simultaneously estimate image content and coil sensitiv¬ 
ities from all data, which minimizes the amount of additional calibration data 
required. Subspace methods do not directly estimate sensitivities, but learn a 
signal subspace from calibration data. These algorithms can adapt to experi¬ 
mental conditions that violate the sensitivity-based signal model formulated in 
Equation [I] For this reason, they are more robust to certain kinds of errors. 

3.1 Non-linear Inverse Reconstruction 

Starting with the signal equation (Eq. [l]), but now considering both image and 
coil sensitivities as unknowns, one obtains a non-linear inverse problem related 
to blind multi-channel deconvolution. Modelling the coil sensitivities as smooth 
functions in a Sobolev space H l (Cl,<C N ), the non-linear version of the forward 
operator can be written as: 


F : L 2 (RC) x H l (Sl,C N ) -> C°°{M. d ,C N ) 

x := (to, Ci, • • • ,c N ) y 


( 11 ) 
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Many auto-calibrating parallel imaging methods reduce the reconstruction prob¬ 
lem to a linear problem by first estimating the sensitivities Cj from a subset of the 
data and then solving for the image using these fixed estimates using a conven¬ 
tional linear reconstruction. Because this is sub-optimal, improved algorithms 
have been developed, which solve the non-linear inverse problem [561 57', I58| . 

In non-linear inversion [53 :, a regularized solution is defined as the solution 
of the minimization problem 


N 

(m, ci, ...,c N ) = argmin x || Fx - y\\ \ + aR(m) + 0^Q(cj) . (12) 

j=i 

Here, a smoothness penalty Q{cj) restricts the solution to smooth coil sensitiv¬ 
ities. The Iterative Regularized Gauss-Newton Method (IRGNM) )5Sj is used 
to iteratively update an estimate of the solution based on a linearization of the 
original problem: 


x n+ i - x n = argmin Ax \\DF Xn Ax + Fx n - y\\\ 

N 

+ a n R(Am +m n ) +/3^^Q(Acj + Cj) (13) 
j =i 

Here, DF Xri is the Frechet derivative of F at the current estimate x n and the reg¬ 
ularization parameters a n ,/3 n are reduced in each iteration step. The smooth¬ 
ness penalty can be chosen as Q(cj) = ||(1 + sA) J Cj||| (with some constants 
s,/). This penalty can be transformed into a / 2 -norm by expressing the sen¬ 
sitivities using Fourier coefficients re-scaled with a positive definite diagonal 
matrix, which avoids bad conditioning of the of the reconstruction problem. 
For quadratic regularization of the image, i.e. R(x) = ||x — xo||l> the algorithm 
then has the explicit update rule 

x n+1 - x n = (DF Xri DF Xn + a n I)~ k (DF Xn (y - Fx n ) + a n (x n - x 0 )) ■ (14) 

Non-linear reconstruction methods can be applied to non-Cartesian sampling mi 
[6T| [62] and extended to include non-linear penalties [321, |63;. 

One limitation of non-linear methods is that they may need an initial guess 
close to solution to converge to the correct global minimum. While it is usually 
sufficient to set the image to a constant value and the coil sensitivities to zero, 
in some cases a guess closer to the true solution is required. In this case, any 
direct estimation method can be used to estimate a set of approximate coil 
sensitivities, which can then be used to initialize the non-linear method. 

3.2 Calibration Matrix 

The calibration matrix is a fundamental tool which can be used to formulate 
many auto-calibration methods. Reconstruction kernels in GRAPPA [53] and 
SPIRiT [M] are null-space vectors of this matrix [65]. The calibration matrix 
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Figure 6: Construction of the calibration matrix: Overlapping blocks of the 
multi-channel k-space become rows of the calibration matrix. 
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Figure 7: Left: Singular-Value-Decomposition (SVD) of the calibration matrix 
reveals signal Vy and nullspace V±. Right: A point-wise eigendecomposition of 
the operator is shown. This operator is derived from the nullspace 

condition V±Rj:f = 0 for each overlapping patch in k-space. The sensitivities 
(here: from an eight-channel coil) appear at each point as a eigenvector to 
eigenvalue one. 
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is a multi-dimensional multi-channel Casorati matrix constructed from fully- 
sampled patches in a calibration area in the center of k-space (see Fig. [b]). It is 
related to the trajectory matrix from singular spectrum analysis (SSA) |66], and 
also to the lag cross-covariance matrix, which can be estimated as £ = jjC h C 
(with M a normalization constant). Because coil sensitivities are very smooth, 
multi-channel signals have correlations in small local k-space patches. This 
implies that the calibration matrix (and the lag cross-covariance matrix) are 
low-rank, i.e. have a small signal space and large null space (Fig. [7]). If the 
calibration matrix is constructed from an incomplete k-space with missing sam¬ 
ples, structured low-rank matrix completion can be used to recover a completed 
matrix, which is the basis of a calibration-less parallel imaging technique known 
as SAKE [571 . 


3.3 ESPIRiT 


Coil-by-coil reconstruction was originally proposed because combination of all 
channels in SMASH-based parallel imaging sometimes caused phase cancella¬ 
tion 168 1 1531 . In combination with auto-calibration coil-by-coil reconstruction 
has a very advantageous side effect: The reconstruction becomes robust against 
certain kinds of inconsistencies - in particular reconstruction in a tight FOV ( 69] . 
The fundamental reason is that the coil-by-coil reconstruction operator does not 
enforce the strict signal model of sensitivity-based reconstruction schemes for¬ 
mulated in Eq. [T] but represents a convex relaxation of this model. 

ESPIRiT is a new reconstruction algorithm which exploits this. Because 
of shift-invariance the null-space condition should be true for all patches in an 
ideal multi-channel k-space /. Let be an operator which extracts a patch 
around a given k-space position k , a least-squares version of this condition is 
then given by 

X>fVlEf%/ = 0. (15) 

k 


This can be further transformed to a convolution-type coil-by-coil operator W, 
which reproduces ideal k-space signals: 


X>F( / - l ii''ii f )%/ = » 

k 

= / 


w 

Wf = f 


(16) 

(17) 

(18) 


Here, M is the size of a single patch. Transforming W into the image domain 
yields an operator which operates point-wise. Because W reproduces 

ideal k-space signals, the image-domain version reproduces the vector of coil im¬ 
ages Cjm at each point r, i.e. [^ 7—1 W^T 7 ] Cj(r)m(r) = Cj(r)m(r). In other words, 
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Figure 8: The power function for Poisson-disc and random sampling for a par¬ 
ticular set of coil sensitivities. Outside of the sampled area and where random 
sampling leaves holes without samples the power function is high indicating 
that recovery of the k-space at this location has a high error when using parallel 
imaging. 


everywhere where the image is non-zero the vector of sensitivities is a point-wise 
eigenvector to the eigenvalue 1 of the operator . The eigenvector and 

eigenvalue maps from a point-wise eigendecomposition are shown in Figure [7] 
Together, these steps form a computational method to extract accurate coil 
sensitivities from the nullspace of the calibration matrix. 

If the k-space is corrupted and does not fit the ideal model, multiple sets 
of sensitivities can appear in affected image regions as multiple eigenvectors to 
eigenvalue one. An extended forward model can take this additional information 
into account. Respective methods offer robustness to certain kinds of errors 
similar to auto-calibrating coil-by-coil methods such as GRAPPA [S5j . 

4 Sampling and Reconstruction in k-Space 

While the formulation of parallel imaging as an inverse problem is a powerful 
conceptual framework, additional theoretical tools are required to understand 
and evaluate different sampling schemes in k-space. For this purpose, a for¬ 
mulation of parallel imaging as approximation in a Reproducing Kernel Hilbert 
Space (RKHS) has recently been developed [21] ■ A RKHS is a Hilbert space of 
functions with continuous (bounded) point-evaluation functionals. This condi¬ 
tion guarantees that sampling is compatible with the norm of the Hilbert space. 
This formulation of parallel imaging yields a unified framework to formulate 
and analyze sampling and reconstruction in k-space. 

By identifying the norm of the signals in k-space with the norm of the corre¬ 
sponding images in L(f2,C), it can be shown that the ideal k-space signals are 
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an RKHS with the matrix-valued kernel 


K st (k, f) = [ drc s (r)c t (r)e- 2 ™-(Z-n . ( 19 ) 

J n 

This kernel captures all local correlations in k-space induced by the sensitivities 
and exploited in parallel imaging algorithms for recovery of missing k-space 
samples. Given this kernel, a standard formula from approximation theory 
can be applied to obtain interpolation coefficients u S j' k for all channels s = 
1 and known k-space samples k £ S for interpolation to arbitrary k- 

space positions: 


N 




( 20 ) 


kes 


With these interpolation coefficients, unknown values in k-space can then be 
recovered from the acquired samples fj(k(ti)) with the interpolation formula 


N 




( 21 ) 


zes 


When no regularization is used in the computation of coefficients the recovered 
ideal k-space corresponds to the best approximate solution defined before. The 
interpolation formulas used in GRAPPA and SPIRIT and similar methods are 
local variants of this formula with an empirical estimate of the ideal kernel 


(Eq. 191. In addition to this interpolation formula, the link to approximation 
theory yields new insights into sampling in k-space. In particular, a point-wise 
error bound in k-space can be derived COJ: 


N 


I/,(£) - /;(£)I 2 < ll/ll 2 • ^(k, £) - £ £ K jt (k, (22) 


les 


P%(k) 


The power function Pj is computed from the kernel K and the interpolation 
functions u and depends only on the coil sensitivities and the sample locations. 
It can be used to analyze the properties of different sampling schemes for parallel 
imaging independent from any actual imaging data. Figure [8] shows the Power 
function for two different sampling schemes computed for a particular set of coil 
sensitivities. 


5 Compressed Sensing Parallel Imaging 

Compressed sensing is based on the idea that randomized under-sampling schemes 
produce incoherent noise-like artifacts in a transform domain which can then 
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be suppressed using denoising to iteratively recover the original signal HUES- 
It exploits the compressibility of the image information, i.e. a sparse represen¬ 
tation in a transform domain, to make the sparse signal coefficients stand out 
from the incoherent noise. Non-linear regularization terms can then be used 
to suppress the incoherent artifacts and recover a sparse representation of the 
image from the under-sampled data. Because MRI acquires data in the Fourier 
domain and is flexible enough to use almost arbitrary sampling schemes, this 
idea can be applied directly mm- 

Parallel imaging can be synergistic-ally combined with compressed sens¬ 
ing [3111331IM1EH- This combination leads to exactly the same optimization 
problems already considered for parallel imaging alone, but requires incoher¬ 
ent sampling schemes suitable for compressed sensing. The most important 
schemes in practical use are variable-density Poisson-disc sampling and radial 
trajectories. Poisson-disc sampling guarantees that samples are not too close 
together. This would waste sampling time, because k-space positions which 
are close are highly correlated and can already be recovered using parallel imag¬ 
ing. Variable-density schemes have several advantages: They equalize the power 
spectrum of the missing samples, provide graceful degradation in case full recov¬ 
ery is not possible, and can be used for auto-calibrating parallel imaging when 
the k-space center is fully sampled. Methods which combine parallel imaging 
and compressed sensing represent the state of the art in image reconstruction, 
as demonstrated by their use in demanding applications such as in pediatric 
imaging without sedation [73 E3 CE2- Figure [9] shows an image from a pe¬ 
diatric patient reconstructed with parallel imaging compressed sensing at an 
acceleration factor of about seven. 


6 Conclusion 

Image reconstruction for parallel imaging can be formulated as an inverse prob¬ 
lem. Based on this formulation, advanced iterative algorithms can be developed 
which (1) make use of optimal (Cartesian or non-Cartesian) sampling schemes, 
and (2) extend parallel imaging with advanced non-linear regularization terms. 
These ideas are combined in recent methods for compressed sensing parallel 
imaging, which currently represent the state of the art in image reconstruction. 
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